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We develop a theory for fluctuations and correlations in a gas evolving under ballistic annihilation 
dynamics. Starting from the hierarchy of equations governing the evolution of microscopic densities 
in phase space, we subsequently restrict to a regime of spatial homogeneity, and obtain explicit 
predictions for the fluctuations and time correlation of the total number of particles, total linear 
momentum and total kinetic energy. Cross-correlations between these quantities are worked out 
as well. These predictions are successfully tested against Molecular Dynamics and Monte-Carlo 
simulations. This provides strong support for the theoretical approach developed, including the 
hydrodynamic treatment of the spectrum of the linearized Boltzmann operator. This article is a 
companion paper to Ref. jl| and makes use of the spectral analysis reported there. 

PACS numbers: 51.10.-|-y,05.20.Dd,82.20.Nk 

I. INTRODUCTION 

Systems where particles may react, change chemical or physical nature and ultimately disappear, model a rich variety 
of phenomena and provide prominent situations to develop and test the foundations of non-equilibrium statistical 
mechanics (see e.g. [1, H, @| and references therein). When reactions are controlled ballistically, the system can 
be modeled by an assembly of hard spheres or disks which annihilate with probability p or collide elastically with 
probability 1 — p everytime two particles meet each other Q . Within the framework of this probabilistic ballistic 
annihilation (PBA) model, most of the work carried out up to now has focused on the kinetic equations for the one- 
body distribution function and the information following from them 0, d, H, [13) El • In particular, the hydrodynamic 
equations, with explicit expressions for the transport coefficients, have been derived by using a generalization of the 
Chapman-Enskog expansion il2j . Our companion paper Jj , where we have established the hydrodynamic description 
in the context of the eigenfunctions and eigenvalues of the linearized Boltzmann collision operator, falls in this vein. 
Within this formalism, the conditions in which the hydrodynamic description is expected to apply are somehow more 
transparent, and can be expressed in terms of some properties of the linearized Boltzmann collision operator. 

In the present paper, the goal is to go beyond the study of one body quantities: the focus is on fluctuations and 
correlations. To this end, we use the tools and ideas developed in the context of the linearized Boltzmann equation. 
The dynamical behavior of the correlations in the dilute limit can indeed be obtained in terms of the linearized 
Boltzmann collision operator. The study of correlations in the PBA model allows to go beyond the description at the 
level of average values, and to characterize how global magnitudes (such as the total number of particles or the total 
energy) fluctuate around their average. It has already been shown in other classes of dissipative systems, such as in 
granular systems, how the knowl edg e of fluctuations is relevant in order to understand the behavior of the system 
when vortices or clusters develop [l3l . [T^ , or even in simpler situations where the system is still homogeneous [15- , 16j . 
The main goal of this paper is to formulate a theory of fluctuations for the PBA model in the dilute limit and to apply 
it to one of the simplest possible state, namely the homogeneous decay state, exploiting its scaling properties. This 
will allow us to obtain explicit expressions for the distributions characterizing the velocity correlations in the system, 
and to compute the statistics of the total number of particles, total momentum and total energy, which decrease 
monotonically due to the annihilation process. 

The paper is organized as follows. In Section |TT1 we present the general framework of the hierarchy method [iTj 
which allows to write the evolution equations of correlation functions. The specific case of the homogeneous decay 
state is considered in Section IIIIl where the scaling properties of this state are used to simplify the equations. After 
briefly recalling in Section IIVI how fluctuations and correlations of global observables can be computed from the 
knowledge of the two-particle correlation functions, we first consider in Section |V] the correlation functions at equal 
time, which give access to the fluctuations of the total number of particles, total momentum and total energy. We 
obtain theoretical predictions for the asymptotic scaling regime as well as for the short time behavior, and we test 
these predictions against numerical simulations (both Molecular Dynamics and Monte Carlo). In Section IVTl we 
generalize our results to the two-time correlation functions and compare also to numerical simulations. Finally, IVIII 
contains some discussions of the results and our conclusions. For the sake of readability, this paper contains some 
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overlap with its companion Repetitions have been kept to a minimum though, and we therefore refer to [1| for 
several technical details. 



II. GENERAL FRAMEWORK 

The system we consider consists of a dilute gas of identical smooth hard spheres or disks of mass m and diameter 
CT, moving ballistically in dimension d. The particles are supposed to undergo only binary collisions. When two 
particles collide they annihilate with probability p or collide elastically with probability I — p- In this probabilistic 
ballistic annihilation (PBA), there is therefore no conserved quantity (except for p — 0) and the number of particles 
decreases steadily. In this section, we will show how to obtain evolution equations for the correlation functions in 
this system. The general idea of the method is to derive a closed set of equations for the distribution functions 
describing the fluctuations by using the same kind of approximations as needed to derive the kinetic equation, in our 
case the Boltzmann equation. In this way, a unified formalism provides the usual kinetic equation as well as evolution 
equations for the one-time and two-time correlations. 

Let Xj = {Rj(t), Vj (i)} denote the position and velocity of particle j in the system at time t. Both, Rj{t) and 
Vj(t) are parametric functions of the initial positions and velocities of all particles. Microscopic one- and two-particle 
densities in the phase space are defined by 

N 

F,{xi,t) = ^5[.Tl-X,(t)], (1) 

1=1 

N N 

F2{xi,X2,t) = ^^(5[xi-X,(i)]5[a:2-X,W], (2) 

1=1 i=£j 

and higher order functions can similarly be defined. Here and in the following, the lower case variables Xi = {r^jWi} 
are field variables referring to a particular point in phase space. 

The initial state of the system is characterized by a point in phase space, T = {Xi, . . . ,Xjv}, which is chosen at 
random according to a probability p(r,0). Introducing the notation (G) = JdrG{T)p{r,0) for the average over the 
initial conditions, the averages of the microscopic densities Fs{xi, . . . ,Xs,t) over p(r, 0) are the usual one-time reduced 
distribution functions 

fs{Xl, ■■■,Xs,t) ^ {Fs{xi, . . ■,Xs,t)). (3) 

Similarly, two-time reduced distribution functions can also be defined in terms of the microscopic densities as 

fr, s{xi, ■ ■ .,Xr,t;x[, . . ■,x'^,t') = {Fr{xi, . . . , Xr , t) Fs{x[, . . .,x'^,t')). (4) 

For simplicity we will consider t > t' > in the following. 

We now introduce the two-particle correlation functions through the usual cluster expansion. The one-time corre- 
lation function 52 and the two-time correlation function hi^i are then defined by 

f2{xi,X2,t) = fl{xi,t)fi{x2,t) + g2{xi,X2,t), (5) 

fi,i{xi,t;x2,t') = fi{xi,t)fi{x2,t') + hi^i{xi,t;x2,t'). (6) 

It is easy to show from the definition of /i, /2 and /i.i that 

hi^i{xi,t;x2,t) = g2{xi,X2,t) + 5{xi - X2)/i(xi,t). (7) 

The case of deterministic annihilation (p = 1) was considered in reference [lO|. The hierarchy of equations for the 
reduced distribution functions is then shown to be similar to the one describing elastic collisions, once the binary 
elastic collision operator is replaced by the operator describing annihilating collisions. In the PBA case, assuming 
molecular chaos, i.e. that no correlations exist between colliding particles^ the equation for Ji{xi,t) is the Boltzmann 
equation 

/i(a:i,t) = J[/i,/i], (8) 
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L^°\x,) = vi- A, (9) 



where 



J[fiJi] ^ y dx2(5(ri2)To(vi,V2)/i(xi,0/i(x2,t), (10) 

and 

ro(vi,V2) = a'^-iy"d^e(vi2-^)(vi2-^)[(l-p)(&-i-l)-p], (11) 

is the PBA binary colhsion operator. In the above expressions Vi2 = Vi — V2 is the relative velocity, ri2 = ri — r2 the 
relative position, Q is the Heaviside function, (t a unit vector joining the centers of two particles at collision, and 
is an operator that replaces all the velocities Vi and V2 appearing to its right by the precoUisional values and 

b^^vi =vl = vi - (vi2 • (t)(t, (12) 

K^^2 = V2 = V2 + (Vi2 • (t)(T. (13) 

The equation for the correlation function g2 can be obtained under the same hypothesis required to derive the 
Boltzmann equation, following the same lines as in reference [isj in the case of inelastically colliding particles, as 



^ + i(°)(xi) + L("Hx2) - K{x,,t) - K{x2,t) 



g2{xi,X2,t) 

= <5(ri2)To(vi, V2)/i(a;i, i)/i(a;2, i), (14) 
where we have introduced the linear operator K{xi, t) 

K{x,,t)^ /da;35(r,3)ro(v„V3)(l+P.3)/i(2:3,i), (15) 



and where the permutation operator Vab interchanges the labels of particles a and b in the quantities on which it acts. 
Finally, the evolution equation for hi^i reads 



^^+L^°\xi)-K{xi,t) 



hi^i{xi,t;x2,t') = 0, t>t', (16) 



that has to be solved with the initial condition ([7|), hi_i{xi,t'; X2, t') = g2(a;i, a;2, t') + 5{xi — X2)fi{xi,t'). 

The equations for the correlation functions /ii.i and g2 contain a part corresponding to free streaming and another 
one which corresponds to collisions. In particular, the evolution of the one-time correlation function g2 due to collisions 
has two parts: one due to collisions of particle 1 or 2 (corresponding to the indices of the correlation function) with a 
third particle, which is governed by the Boltzmann collision operator; and a second one, due to collisions of particle 
1 with particle 2, which can be written in terms of the one particle distribution function as a consequence of the 
molecular chaos hypothesis. In fact, as in the case of the inelastic granular gas [T^, the only difference between the 
evolution equations of the correlation functions for the PBA and for a system of elastic particles lies in the substitution 
of the elastic binary collision operator by the operator for the PBA model, 10(^1, V2). However this does not give 
any a priori guarantee on the range of validity of these equations. This prescription, i.e. how small the density of the 
system must be so that the above kinetic equations provide an accurate description, might depend on the parameter 
p, and also on the particular state being considered. 

We will see in Section IIVI how the knowledge of and 32 allows to obtain the correlation functions of any 
observable which is a function of the particles positions in phase space. 

III. HOMOGENEOUS DECAY STATE 

As recalled in the previous companion paper [l| , the Boltzmann equation for the PBA model ([5]) admits a particular 
solution //f(v, i) describing a spatially homogeneous decay state (HDS), in which all the time dependence is contained 
in the evolution of the density nnit) and the temperature T/f (i), which are defined as in standard Kinetic Theory as 



nnit)^ JdyrfH{yr,t), ^njf(i)r^(t) = j dy^'^v^ fni^^t). (17) 
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Although there exists no rigorous proof of its stabihty nor of the fact that such a state should be approached from 
arbitrary initial conditions, numerical results obtained by Molecular Dynamic simulations and by the direct Monte 
Carlo method strongly support the existence of such a homogeneous solution [j, 9]. In this section, we review for 
completeness the evolution equation of the one-particle distribution function and obtain the equations for adequately 
rescaled correlation functions. All quantities concerning this homogeneous decay state will be labeled by an index H. 
In the HDS, the one body distribution function does not depend on space and follows the scaling form [T3| 



(18) 



where nn (t) is the uniform density, vh (t) is the thermal (root-mean-square) velocity related to the granular temper- 
ature Tnit) by 



VHit) = 



2THit) 



1/2 



(19) 



and x(c) is an isotropic function depending only on the modulus c = |c| of the rescaled velocity c ~ v/u/f(t). 
Moreover, the density and temperature fields evolve according to (T^ 



dt 
dTnit) 
dt 



(20) 

(21) 



where VH{i) is the collision frequency of the corresponding hard sphere fluid in equilibrium (with same temperature 
and density) 



yH(t) 



nH{t)Tll\t) a 



ml/2 (d 2)r(d/2) ' 

and the dimensionless decay rates C„ and Ct are functionals of the distribution function: 

yiici Jdc2T{ci,C2)xH{ci)xH{c2), 
-2cl 



1 

2 



- 1 T(ci,C2)Xff(ci)XH(c2). 



(22) 

(23) 
(24) 



In these expressions, 7, which does not depend on time, is given by 7 = 2nH{t)vH{t)cr'^ ^/vH{i) = {d 
2)v^r(d/2)/ (47r(''"i)/^), and the binary collision operator T(ci,C2) takes the form 



T(ci, C2) = Jd&e{ci2 ■ o-)(ci2 • 6-)[(l - p)b-^ - 1] 



Finally, the scaled distribution function xh{c) obeys the equation 

d 



P 



(dCr - 2C„) + Ctci 



dci 



Xh(ci) = 7 jdC2T{ci,C2)XH{ci)XH{c2)- 



(25) 



(26) 



The operator h'^^ in the last equation is defined again by equations and p^ . but substituting (vi, V2) by (ci , C2). 
The analytical form of xh is not known, but its behavior at large and small velocities has been determined [9|, llOj ] . 
As in the companion paper [l|, we will use here the approximate form of the distribution function in the so-called 
first Sonine approximation, which is valid for velocities in the thermal region, and all the functionals of X-f/(c), like 
the decay rates and the transport coefficients, will be evaluated in this approximation [7|, |9|. 

Considering the scaling form for the one-particle distribution function, it is convenient to introduce the rescaled 
correlation function gn through 



g2,H(Xl,X2,t) = 2dU\ 9H\T, ri2,Ci,C2), 



(27) 
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where we have taken into account that the system is invariant under space translation, so that 52,// depends on 
ri2 = Ti — r2 and not on ri and r2 separately. The dimensionless time scale r 



1 

^=2 



t 

dt'i^nit'), (28) 



is proportional to the number of collisions in the time interval [0,t\. The equation for the reduced function gn in 
these units reads then 



+ A(ci) + A(c2) - 2K« - Ih{t)ci2 



gniT, ri2,Ci,C2) 



or or 12 

- -5(ri2)7T(ci, C2)xh(ci)xh(c2), (29) 

where we have also introduced the length scale I nit) = 2v}j{t) / V}j{t)^ which is proportional to the instantaneous 
mean free path, and the linearized Boltzmann operator (see previous paper) 

A(c,)/i(c,) = 7 JdC3T{c„C3){l+V^3)XHic3)h{c^) 

+ p(2C„ - dCrMc,) - pCtc. • -^hic). (30) 

OCi 

Similarly, we define a rescaled two-time correlation function hn through 

hi^i.H{xi,t;x2,t') = -j-^^;^^^/iH(ri2;ci,r;c2,r'), (31) 
and obtain the following evolution equation: 

d 



^ftif(ri2;ci,r; C2,r') 



A(ci) -;H(r)ci 

OTi 



/iff(ri2;ci,T;c2,r'), (32) 



with hH{ri2]Ci,T\ C2, r) = ^^(t, ri2, Ci, C2) + 5{ci - C2)(5(ri2)xH(ci). 

It is interesting to note how, in this representation, all the time dependence due to the reference state is absorbed 
in the free streaming term through the function Ih{t), proportional to the mean free path. The evolution of the 
correlation functions moreover will be determined by the properties of the linearized Boltzmann operator A, which 
we have already studied in the companion paper [l| and which we will recall in Section [V] 

IV. FROM PARTICLE CORRELATION FUNCTIONS TO CORRELATIONS AND FLUCTUATIONS OF 

GLOBAL MAGNITUDES 

In this section, we will show how the previously presented framework for correlation functions will allow us to study 
the fluctuations and correlations of global quantities for a PBA system in the homogeneous decay state. In particular, 
we will focus on the total number of particles N, the total momentum P, and the total energy E. 

Consider indeed two dynamical variables of the form 

N 



^[r(t)] =^a(V,) = fdxiaiv,)Fi{xi,t) 

N 

B[m] = I^KVO = jdx2h{w2)Fi{x2,t) (33) 

2 — 1 

where a and h are functions of the particles' velocities V^, and Fi{xi^ t) is the microscopic density in phase space ([ij. 
Taking a = 1, a = v and a — mv|/2 yield for A the total number of particles N, the total momentum P, and the 
total kinetic energy E, respectively. The deviations 6A{t) = A{t) ~ {A{t))H and SB{t) = B{t) - {B{t))H of A or B 
from their average values in the HDS (denoted by (...)//), define the fluctuations of both magnitudes, which have 
average zero and correlations 

{SA{t)SB{t'))H = {A{t)B{t'))H - {A{t))H{B{t'))H. (34) 
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It is then straightforward to use the definition of the two-time correlation function /ii.i in Eq. ([6]), to obtain 

{SAit)5B{t'))H = jdTtjdwtjdY2jd^2a{^i)b{^2)hi,i,H{vu-fut;r2,^2.t'). (35) 

In particular, for t = t', this leads to 

{SA{t)5B{t))H = V y"dva(v)&(v)/H(v,t) 

+ V jdvijd^f2a{wi)b{w2) yrfri2g2,ff (ri2, Vi, V2, t), (36) 

where V = J dri is the total volume of the system. These formulas show how the correlations of two different global 
magnitudes are determined by the one particle distribution function and by the correlation functions. The one particle 
distribution function is known in the HDS in the first Sonine approximation 0, 

V. FLUCTUATIONS IN THE HDS 

Let us focus in this section on the one-time correlation function gn- We will only need functions a and b which depend 
on the velocity degrees of freedom, so that it is convenient to integrate out the spatial dependence by introducing 

(/)h(t,Ci,C2) = Jdri2gH{T,ri2,Ci,C2), (37) 



whose evolution is obtained from (1291) as 

A_h-(t, ci,C2) = -jT{ci,C2)xh{ci)xh{c2)- (38) 



d 

-— + A(ci)+A(c2)-2K« 



Given an initial condition (0, Ci, C2), this equation (|38p can be formally integrated as 
0ff(T,ci,C2) = e(^('=i)+^('=^)-2K„)-^^(0,ci,c2) 

Jo 

= e(A('=^)+^('=^)-2K.K[^^(O,Ci,C2)-0|,(Ci,C2)]-f <^|,(Ci,C2) 



(39) 



where (/)Jj(ci,C2) is the solution of 

[A(ci) + A(C2) -2K«]0|/(ci,C2) = -7r(ci,C2)xH(ci)x//(c2), (40) 

where we implicitly assumed that A is invertible. This happens to be the case, see below. The spectral properties of 
the linearized Boltzmann operator A are thus crucial for the evaluation of (j)H ■ We therefore start by recalling these 
properties. 

A. Spectral properties of A 

In our companion paper we have analyzed the eigenvalue problem associated with the linearized Boltzmann 
operator A 

A(c)e^(c) = Xp^fsic). (41) 

We have in fact restricted ourselves to the hydrodynamic part of A, defined, by those eigenvalues of the balance 
equations for the number densitVjmomcntum, and temperature following from the homogeneous linearized Boltzmann 
equation. Such eigenvalues are [l2'| 

Ai=0, A2 = -p(CT + 2Cn), A3 =Kt. (42) 
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Although we were not able to prove that these eigenvalues are indeed the hydrodynamic ones, the self-consistency 
of the resulting description and the successful comparison with numerical simulations have validated this assumption 
[1] . In the previous paper we also obtained the corresponding eigenfunctions 

6(c) = Xff(ci) + ^ • [cxh(c)], (43) 

Uc) = zxh{c) - ^ ■ [cxh(c)], (44) 
oc 

6(c) ^ (45) 
oc 

with z = a function of the probability of annihilation p. The eigenvalue A3 is c?-fold degenerated and we denote 
^3i, i = 1, d the corresponding eigenvectors. The scalar product of two functions /(c) and g{c) is defined as 

(f\g)^ JdcxH\c)r{c)g{c), (46) 

/* being the complex conjugate of /. The eigenfunctions ^/j given in (P51) - (|l5)) are not orthogonal, as a consequence of 
the operator A being non-Hermitian in the associated Hilbert space. On the other hand, it is easily verified that the 

set of functions {Ci; 6; I3} = |xh(c) - (5 + ^) X/f (c); (^i + 4) Xff (c); cxh(c)| verify the biorthogonality 

condition (^/3|^/3') = Sp.p', for P, j3' = 1,2,3. 

The eigenfunctions of the operator [A(ci) + A(c2) — 2p(„] that appear in equation PO)) are then easily seen to be 
C/3i(ci)C/32(c2), with 

[A(ci) + A(c2)-2K„]?/3i(ci)e&(c2) = (A^, +A^, -2K„)C/3i(ci)^&(c2). (47) 

Since Cn > Ct p^ . and under the assumption that the norm of the "non hydrodynamic" eigenvalues are always 
greater than the hydrodynamic ones, the eigenvalues of [A(ci) + A(c2) — 2pC„] are therefore all negatives. This 
has the important consequence that the exponential term in (|39p decays to zero and that the large time limit of 
(/)_f/(T, Ci, C2) is 01^(01,02), solution of Eq. ([iO)) . 



B. Hydrodynamic part of the correlation functions 

Obtaining the full spectrum of A is a formidable task. We will here assume, as in the companion paper [T], that 
the kinetic (non-hydrodynamic) modes have a fast decay, and work in the hydrodynamic subspace spanned by the 
functions S,fj defined in the previous subsection. In that purpose, we generalize the definition of the scalar product 
given in (|46p to two-velocity functions by 

(/(C1,C2)|.9(C1,C2)) = JdCiJdC2XH^ici)XH^{C2)f*{Cl,C2)gici,C2). (48) 

This allows to define a projector operator P12 onto the space spanned by the functions '?/3i (ci)^;32 (02) as 

3 3 

Pl2/(C1,C2) ^ E E (eft(ci)Cfe(c2)|/(ci,C2))Cft(ci)efe(C2), (49) 

/3i = l/32 = l 

and the "hydrodynamic part" of 0//(r, Ci, C2) and 0|f(ci,C2) are by definition 

3 

(?!)^'(t,Ci,C2) = Pi2(/'h(t,Ci,C2) ^ A/Si ,fe (■^)?/3i (ci )^/32 (C2 ) , (50) 

/3l,/32 = l 
3 

'/>H^Hci,C2) = Pi2'?!)h(Ci,C2) = a^i,/32^/3l(Cl)C/32(c2)- (51) 

/3i,/32 = l 

We can now obtain a closed equation for 0^'' by applying the operator P12 on both sides of equations ([39|) and 
(HOI) , under the additional assumption that 

Pi2A(c,) = Fi2A(c,)Fi2. (52) 
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A theoretical estimation "a priori" of the accuracy of this approximation would require to know more than it is 
available at present about the spectrum of A(c) and its adjoint. Therefore, it will be taken as a working hypothesis 
to be validated later on by comparing the predictions it leads to with the results from numerical simulations of the 
system. It is worth emphasizing that since A leaves the hydrodynamic subspace invariant, Equation (|52p is equivalent 
to the commutation relation P12A = AP12. Proceeding further, we obtain from pQ]) 



(r, ci, 02) - e(^(^i)+^(^^)-2K„)rp^2 [^^(0, d, 02) - 0?^ (ci, C2)] + Pi2<^?^ (ci, C2) 

Pi, 02 = 1 



(53) 



where we have introduced Ap-^ = ap^ (0) — a^^ . We show in Appendix 1X1 how to obtain explicit formulas for the 
coefficients a^^ in terms of the cooling rates, ('n and Ct, and other coefhcients which are also functionals of the one 
time distribution function xh- The values of ap^^p^iO) depend on the initial condition i/i/f (0, Ci, C2). For the specific 
case in which the variables TV, P and E do not fluctuate at t = 0, and taking into account that the system is in the 
HDS, the coefficients ap^^p2{0) are calculated in Appendix [Bl 



C. Hydrodynamic approximation for global fluctuations 



In this section we compute the correlation functions of the global observables by replacing (j)H implicitly appearing 
in (|36p by its hydrodynamic part, (/i^^. This can be done invoking the relation 



(^"A(ciKft(c2)|/(ci,C2)) = (^;3,(ci)Cft(c2)|/('')(ci,C2)), 



(54) 



for Pi = 1, 2, 3. However, it must be stressed that the theoretical prediction for 0^'' in Eq. ((551) . has been calculated 
using the approximation (|52p . 

If we substitute a(v) = 1 and 6(v) = 1 in (1551) . we obtain for the fluctuations of the number of particles 



{6N^t))h = Nh{t) 



Jdcxnic) + Jdci y"dc2(/)^'' (r, €1,02) 



(55) 



where we have introduced the notation Nh = {N)h- In order to calculate the fluctuations of the total momentum we 
substitute a(v) = Vi and fe(v) — Vj in equation (|36p and obtain 



JdcciCjXnic) + y"dciydc2CiiC2j0^^(T, Ci,C2) 



For the energy we substitute a(v) = 6(v) = ^mv^ so that we have 



{5E\t))h^—Nh{t)v%{t) 



jdcc^XHic) + jdCi JdC2clcl(j)^^\T,Ci,C2] 



(56) 



(57) 



Finally, we can calculate the correlation between SN and SE by taking a(v) — 1 and b{v) 



{SN{t)SE{t))h = -^NHir)v%iT) 



Jdcc^XH{c) + Jdci Jdc2cl(j)'-^\T, 01,02) 



After some algebra, we obtain 



(58) 



{6N^t))h = Nh{t) \i + + 2zaJ 3 + z^a^^ + ^lae-^PC-- + 2zAi,2e-P(^^+4^")" + z^A2,2e-^''^^^+''^"^^] (59) 



{SP{T)dP,{T))H = d.,NH{T)vUr) \ 2 +a3.,3. 



1 - e 



(60) 



{SE^r))H^—NHirKir) 



did + 2 



^ (1 + a2) + -al, d\l + -)al2 + + + -A.^.e-^^^-^ 



-d' (1 + A,,2e-P«-+^'^-^^ + d'[l + A2,2e 



-2p(Ct+3C„)t 



(61) 
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and 



{5NiT)SE{T))H = -TrNHiT)vUr) 



d d , 

+ dAi 26 



dal 



-p(Cr+4C„)r 



(l + I) "2,2 - 
dz + ^2^2 



:^i,ie 



-2pC«T 



-2p(Ct+3C„)t 



(62) 



where 02 is related to the fourth moment of xh{c) through J dcc'^XH{c) = — ^[1 + a2], and has been evaluated in 
the first Sonine approximation in p^ . All the functions ^ and are evaluated in the Appendices lAl and [b1 

At this point, it is important to note that the equations l|59 p - (|62p have been obtained under the assumption that 
the system is in the homogeneous decay state at all times, i.e. that the one particle distribution function is Xiiic) 
for all the time evolution. For p < 1, if we start with an arbitrary initial condition, numerical simulations show that, 
after a few collisions, the distribution function reaches the scaling regime given by equation IjlSp. Then, the evolution 

of (^^^ is given by ([55]) and one expects that the same correlation functions ([5^ - ([S^ will be obtained in the long time 
limit, independently of the initial condition. This will be confirmed in the next section by numerical simulations. 

Equations (|59p -(|62 p lead to a certain number of theoretical predictions. In particular, they imply that the ratios 
{5N\t))/Nh{t), {6E^t))/{Nh{t)v%{t)), {5Pf{T))l{NH{T)vl{T)) and {5N{t)5E{t)) /{Nh{t)vI{t)) reach station- 
ary values at large times. The approach to these stationary values is exponential in r, and is slower for the correlations 
of the total momentum, since the argument of the exponential is p{Cn — Ct)t^ while the other quantities evolve on 
faster time scales. 



D. Numerical simulations 



We now compare our theoretical predictions with the results of Molecular Dynamics (MD) and Direct Simulation 
Monte Carlo (DSMC) of a freely evolving system of N hard disks of diameter a which annihilate with probability p 
or collide elastically with probability 1 — p everytime two particles meet each other. In the MD case, the particles 
were localized in a square box of size L with periodic boundary conditions. The event driven algorithm [l^ has been 
used and the initial density has been chosen low enough to be always in the dilute limit. The parameters for all the 
MD simulations were N{0) = 10^ nij(0)cr2 = 0.05, Ti/(0) = 1 and < p < 1. In the case of the DSMC simulations 
we have used Bird's algorithm [13] with the same values of the parameters, except the density that plays no role. The 
initial velocity distribution is a Maxwellian in both cases. We have measured the time evolution of the total number 
of particles and the total energy, averaging the data over various initial conditions (the total momentum fluctuates 
around zero). We have first checked that the equations (|20p - (PT1) . with the theoretical predictions derived in for 
the cooling rates, correctly describe the decay of the average global quantities. In the same way, we have obtained 
the averaged values of N'^{t), E'^{t), Pi{t) and N{t)E{t) (the correlations between Pi and N 01 E are zero). 

Figures [U [2] and [3] show the time evolution of the various one-time correlation functions considered, for p = 0.5 and 
p — 0.8. The DSMC results have been averaged over 4000 trajectories while the MD simulations has been averaged 
over 150 trajectories (the DSMC method being computationally less expensive, it is then possible to average over a 
larger number of initial conditions than for the MD simulations). The dashed lines are the theoretical predictions, 
equations (|59p - ([5^ . Note, however, that the system is not initially in the HDS : the initial distribution function is 
a Maxwellian and not xh- Nevertheless, as the difference between these two distributions is very small (at least for 
thermal velocities since 02 0.1) and as the stationary values depend very weakly on p, equations (|59p - (|62p predict 
quite well the time evolution measured in the simulations. The ratios {SN'^{t))/Nh{t), {5E'^{t)) / {Nh{t)v^{t)), and 
{5N{t)5E{t)) /{NH{T)v'jj{T)) reach stationary values as predicted. The fluctuations of the total momentum evolve 
more slowly, as also predicted, and the stationary value of {6P^ (t)) / {Nh {T)v'j^ (t)) is barely reached. Note that r = 4 
corresponds for p = 0.5 to a total number of particles at the end of the simulation N ~ 1700. 

We have performed simulations starting with other initial conditions further from the HDS. The initial velocity 
distribution function has been chosen as a constant function in a square centered in the origin in the velocity space 
such that the initial temperature is unity. As seen in Figures |4] and [5l we obtain a different short time evolution but 
the scaled moments still converge towards the HDS values, that are represented by the dashed lines. The convergence 
is slower as we increase the value of p, and {5E'^{t)) / {Nh{t)v'^{t)), the magnitude which depends on the higher 
moments of the velocity distribution, is the most affected. 

Figures [6] and [7] show the comparison between the stationary values of the various ratios measured in the simulations 
and the theoretical predictions in Eq. (|59ll62p at large r. The agreement is very good for all values investigated. We 
have also computed the probability distribution for the number of particles, energy and momentum. As we can see 
in Fig. [51 where we have considered a system with p = 0.5, they are correctly described by a Gaussian distribution. 
The figure displays the distribution at four different times, showing that the shape of the probability distributions 
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Figure 1: Scaled second moment of the fluctuations as a function of r for a system with p — 0.8. These results are from DSMC 
simulations and have been averaged over 4000 trajectories. 
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Figure 2: Second moment of the fluctuations of the number of particles (left panel) and of the total energy (right panel) as 
a function of the number of collisions per particle r, for a system with p — 0.5 and initial number of particles = 10^. The 
dashed lines are the theoretical predictions. 



does not vary during the dynamical evolution. Similar results have been obtained for the probability distribution of 
the total momentum. 



VI. TWO-TIME CORRELATION FUNCTION IN THE HDS 



In this section, we study the two-time correlation function of the global quantities in the HDS. To this aim, we 
consider two dynamical variables A{t) and B{t) as in and compute the correlations {5A{t)5B{t'))H for t > t' , 
which are obtained from /ii.i,// through equation (j35p . 

As in the previous section, we start by integrating out the spatial degrees of freedom and consider 

iPh{t, t' ,Ci,C2) = y dri2/iH(ri2;ci,T;c2,r'), (63) 

whose evolution equation is obtained by integrating (|32p over space variables 

d 

— V'h(t,t',Ci,C2) = A(ci)V'if(T, t',Ci,C2). (64) 

This equation has to be solved with the initial condition 

4>h{t' ,t' ,Ci,C2) = Xh(ci)(5(ci - C2) + (/)i/(T',Ci,C2), (65) 
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Figure 3: Correlation between the fluctuations of the total energy and total number of particles (left panel), and second moment 
of the fluctuations of the y component of the total momentum (right panel) , as a function of the number of collisions per particle 
r, for a system with p = 0.5 and initial number of particles A'' = 10^. The dashed lines are the theoretical predictions. 
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Figure 4: DSMC results for the scaled second moment of the fluctuations as a function of t for systems with p — 0.2 (left 
panel) and p = 0.5 (right panel). The dashes lines are the theoretical predictions for the stationary values. The initial velocity 
distribution at r = is uniform in a square domain. 

where we have taken into account the scaling of fn (|18p and g2,H (|27|) . Then, using the approximation ()52[) . we obtain 

Pi2^/'h(t,t',Ci,C2) = e^('=i)(^-^')Pi2^H(T',r',Ci,C2) 

= (6 (ci ) (r , r' , ci , C2 )) Ci (ci ) 

+ {UciMnir', r', ci, C2))6(ci)e-^('^-+2';")(^-^') 

+ Y.(^3^{ciMHir',T',c,,C2m3,{ci)eP'^^'-^-^'\ (66) 

i 

In the large time limit r, r' (X), r — r' finite (and positive) we can replace i1^h{t\t\ci^C2) by Xh(ci)<^(ci — C2) + 
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Figure 5: Same as in Figure |4] but for a system with p = 0.1 
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Figure 6: Average steady state of the scaled second moment of the fluctuations of the number of particles and of the energy 
fluctuations as a function of the annihilation probability p. The dashed lines show the large r predictions of Eqs. ((59} and (|61|l . 



6|^(ci, C2), so that 



oPCtC-t-t') 



dm 



(67) 
(68) 

(69) 
(70) 

In the T scale, it can be seen from equations ((20)) and ((2T|) that A'^/f and u// decay exponentially. For A,B = N, E, P, 
the normalized correlation functions 



[-{Z + 2)2^2,2 + {Z + 2)Ai,2] e-^'C^T+SCOlr-r')! 

(<5iV(T)<5ii;(r'))if = TO^7VH(r)«|,(r'){Ai,i-(z + 2)Ai,2 

+ hz(z + 2)^2,2 + zAi,2] e-f J . 



Cab{t,t') = {5A{t)5B{t'))I{5A{t')5B{t')), 



(71) 



become therefore time-translation invariant, i.e. functions of r — r' once the stationary regime for the ratios such 
as {SN'^{t)) /Nh{t) has been reached (see previous section). We have checked numerically that this is indeed the 
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the probability of annihilation p. The dashed line shows the large r prediction of Eq. (|62p . 
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Figure 8: Normalized distribution of the relative energy (left panel) and number of particles fluctuations (right panel) for a 
system with p = 0.5. The symbols are from DSMC simulations and for four different values of r, t\ — 0.69, T2 = 0.98, T3 = 1.37 
and T4 = 2.05. The solid line is a Gaussian with unit variance. 

case, and we compare in Fig. l9l and \T0\ the evolution of Cnn{t — r'), Cee{t — r') and Cne{t — r') measured 
in DSMC simulations (for p = 1, averaged over 4000 trajectories) with the theoretical predictions. The agreement 
is very good. Figure [TT] also shows the theoretical prediction for the decay of the momentum correlation function 
Cpp{t,t') = {6Pi{t)SPj{t')) / {SmT')SPj{T')) for p = 0.5. The characteristic decay time of Cpp is of the order of 
r ~ 4. Because the time to reach the stationary regime for {5Pf {t)) / {N H{T)v'jj{T)) is t ~ 4, as shown in Fig. [31 we 
would need to reach r ~ 8 in the numerical simulations in order to display numerical results for Cpp, which means 
that we would need to consider simulations with an initial number of particles of the order of iV(0) ^ 10^. 
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Figure 9: Decay of the two-time correlation of the number of particles, Cjvjv, and the energy, Cee, for a system with p = 1, 
measured with DSMC simulations. The dashed line is the theoretical prediction. 
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Figure 10: Decay of the two-time correlation of the number of particles and the energy, Cne for a system with p — 1. The 
dashed line is the theoretical prediction. 




T-X' 

Figure 11: Theoretical prediction for the decay of the momentum correlation function Cpp{t, r'), defined in the main text, as 
a function of r — r' for a system with p = 0.5. 

VII. CONCLUSIONS 

In this paper, we have formulated a general theory for fluctuations and correlations in a dilute probabilistic ballistic 
annihilation system. The theory has been particularized to the homogeneous decay state, taking advantage of its 
scaling properties. For this state we have focused on the study of the fluctuations of the total number of particles, 
total momentum and total energy, evaluating the two-time correlation functions between these quantities in the 
hydrodynamic approximation. The fluctuations of the total number of particles, total momentum and total energy, 
once conveniently rescaled, converge to stationary values. The convergence is exponential in the natural time-scale 
T given by the number of collisions, and the corresponding rates are simple combinations of the cooling rates. The 
stationary values are obtained as functionals of the distribution function and can be computed in the first Sonine 
approximation. We have also obtained theoretical expressions for the two-time correlations of global observables. 
All our theoretical predictions have been successfully compared with the results of Molecular Dynamics and DSMC 
numerical simulations, providing strong support for the theory developed here, including the hydrodynamic description 
in terms of the lowest order eigenfunctions and eigenvalues of the linearized Boltzmann collision operator. 

As a side-remark, we note that the correlation functions contain two parts: one coming from the one-particle 
distribution function, and another one that comes from velocity correlations. Nevertheless, it must be stressed that 
the existence of these velocity correlations does not imply a violation of the "molecular chaos" assumption that 
underlies the Boltzmann equation. This is because the latter only refers to the precoUisional part of the two-body 
distribution function (at contact). 

The fact that the two-time correlation functions decay on a time-scale determined by the cooling rates reflects 
the intuitive notion that their dynamic is essentially of macroscopic character, compatible with Onsager's regression 
hypothesis (see e.g. [HI). To analyze this point in a deeper way, we show in Appendix [Cl that a description of the 
fluctuations SN, 6E and 6P in terms of linear Langevin equations can be obtained, using for the deterministic part the 
evolution equations for a linear perturbation around the HDS. The (Gaussian and delta correlated in time) noise terms 
in the Langevin equations can then be adjusted in order to obtain the same amplitudes for the one-time correlation 
functions as with our theory. In this respect, the results derived in Appendix [Cl mav be envisioned as formulating 
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a fluctuation-dissipation theorem for the homogeneous decay state under scrutiny in this paper. The amphtudes of 
the noise terms are however comphcated functions of moments of the one-particle distribution functions, and arc not 
clearly related to macroscopic quantities such as the cooling rates. 
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Appendix A: EXPRESSIONS FOR 



(32 



In this Appendix we compute the expressions for the coefficients a^^ Applying the projector P12 to 
under the assumption Pi2A(ci) = Pi2A(ci)Pi2, 

[A(Ci ) + A(C2 ) - 2Kn] 4>T^ (Cl , C2 ) - -7i^l2T(Ci ,C2)XH (Cl )XH (c2 ) ■ 

Using the definition (|5ip of the coefficients a^^ , we then obtain the set of equations 
3 3 

mi^ft>/92(V + - 2Kn)C/3i(cl)C/32(c2) = -Pl27T(ci , C2 )xH (Cl )xh (C2 ) , 

Pi f32 

and it is straightforward to write 

^ (C/3i(cl)C/32(c2)|7T(ci,C2)xg(ci)xg(c2)) 

Given the expression of the functions {fi(c)}f^]^, and taking into account the relations 

PCn 



we obtain 



2 y"'=^Ciy"'^'^2T(ci,C2)x//(ci)xH(c2), 
I Jdc,Jdc2 (^^-1^T(ci,C2)xh(ci)xh(c2), 



1 2 



1 - 



2(1 + z) 



z J 1 + z 



1 - 



2(1 + z) 



Hp) 



1 



pCt + 4pC 

+P{Cn + Ct) 



2(1 + ^) 

z 



- 1 



2pCn (1 + Z)^ ■ 

p{2Cn + Ct) 



(1 + ^) 



Hp): 



1 



2p(Ct + 3Cn) 



2(l-hz)2 ^^^(l+z)2 

p(3C„ + 2Ct) , bip) 



2(1 + ^)2 (1 + ^)2 

2(Ct ~ Cn) 



where 



Hp) = 7 dCi dC2-^T{Ci,C2)XH{ci)XHic2), 

Hp) = 7 jdcijdc2XH{ci)xH{c2)jda-Q{ci2 ■ a-){ci2 ■ a-)cixC2x- 



|40)) yields, 
(Al) 

(A2) 

(A3) 

(A4) 
(A5) 



(A6) 



(A7) 

(A8) 
(A9) 

(AlO) 
(All) 
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These two functions have been evaluated in first Sonine order with the result 

16(-1 + 4d{d + l))p + a2[256 - 255p + 4d(-64 + (71 + 7d)p)] ^(d_i)/2 



(-16 + 5a2) M_iW2 
32\/2dr(d/2) 



(A12) 
(A13) 



Appendix B: EXPRESSIONS FOR a^Si ^(O) 

In this Appendix we evaluate the coefficient a/3i,/32(0) for the specific case in which we have 

(JAr2(o))=0, (,5P,(0),5P,(0)) =0, 
((5i;2(0)) = 0, ((57V(0)(5i;(0)) = 0. 

Taking these relations into account, it is straightforward to obtain 

-1, 

2 ' 
d 

~2' 
rf(<i+2) 



dci y"dc2 0ff(O, Ci,C2) 

dCi y dC2CiiC2j4'H{0, Cl, C2) 
dCiJdC2cl(j)HiO,Ci,C2) 
dCi JdC2clcl(j)H{0,Ci,C2) 



(l + a2). 



With these expressions and the definition of we get 

1 



ai,i(0) 

«1,2(0) 

02,2(0) 
031,31(0) 



(1 + Z)2 
1 

2(l + z)2 

1 

"(1 + Z)2 
1 

^2' 



-(z + 2) - i(z + 2)2 - ^z2(^ ) 
2^ ^ 4^ ^ 4d ^ ^ 2; 



(2+2)+^<l + 02) 



3 d + 2^^ 

4 + ^(1 + "^) 



(Bl) 
(B2) 



(B3) 
(B4) 
(B5) 
(B6) 

(B7) 
(B8) 
(B9) 
(BIO) 



Appendix C: LANGEVIN EQUATIONS FOR THE GLOBAL MAGNITUDES 

In this Appendix, we will show that it is possible to find a Langevin description for the fluctuations of the global 
magnitudes of the system. The idea is to assume that the global magnitudes obey some equations that can be 
decomposed in a "deterministic part" , which is identified with the macroscopic equations for a linear perturbation of 
the HDS, plus a Gaussian white noise. Because of formulas ([551) - (I5^ let us study the equations for the magnitudes 



SNir) 



5N{t) 



5E{t) = 



5P^{T) 



N]l\T)vH{r) 



A5E{t) 



dmN]l'^{T)vjj{i 



(Cl) 



in order to deal with processes with time independent variances. 

Let us start with the easiest one, the equation for SPi{T). We can define the function 



Wi,k=o(T) 



1 



nH{t)vH{t) 



drjdvvi6f{r,v,t), 



(C2) 
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where 5f{r,v,t) = /(r,v,t) — fH{v,t), with fHiVjt) the distribution function in the HDS. Then, if the generic 
distribution function /(r, v, t) is close enough to the HDS one, the hnear equation for LOi ]s^^o(T) is (see the companion 
paper) 



d \ 

^ - Kt I UJiM=o{T) = 0. 

Then it is straightforward to see that the equation for the macroscopic deviation SP^^ would be 

d 



di 



+ P{Cn~CT) 



(C3) 



(C4) 



where the superscript "M" denotes macroscopic. Now, let us suppose that the equation for the fluctuating SP is of 
the form 



d_ 



+ P(Cn-CT) 



SP^ir) = Rpir), 



with Rpir) a Gaussian white noise with the following properties 

{Rp{t))h - 0, {RpiT)Rp{T'))H = TpSir - t'). 



(C5) 



(C6) 



That is, we consider that the equation describing the dynamics of the fluctuations can be obtained from the macro- 
scopic equation describing the evolution of the system. Under these hypothesis we can calculate {SP{t)SP{t')). The 
solution of equation (jCSP in the long time limit is 

SP,{t)= r dT'e-P^<"~^^'>^^"^'^Rp{T'), (C7) 
Jo 

and the autocorrelation function is, for t > t' 



{5P{r)SP{r')) 



P g-p(Cn-CT)(T-T') 



1 -e 



MCn - Ct) 

We are interested in the limit, r' — > oo, t — > oo, t' — t ^ finite. In this limit we obtain 



{SP{t)SP{t'))h 



MCn - Ct) 



-p(Cn-CT)(T-T') 



(C8) 



(C9) 



Now one can relate this result to the one obtained in the previous section, equation (|68p . that can be expressed in 
our variables as 



{SP{r)SP{T')}H^ l+a 
Comparing equations (|C9p and (jClOp it is seen that if 

Tp = 2p(C„ - Ct) 



3i,3i 



-p{<,^-CT)iT-r') 



-+a. 



3i,3i 



(CIO) 



(CU) 



the Langevin equation (jC5[) is in agreement with the results obtained in the previous section. 

Now we will sketch the derivation of the Langevin equations for the other fluctuating quantities. First of all, we 
are going to start from the macroscopic equation for po and eo defined as 



^0 = ^1^=0(1-) = [dr [dvSf{r,v,t), 

nH[t) J J 



£0 = ek=o(T) 



Jdrjdv^mv^5f{r,v,t). 



These equations are 



dr 



dnH{t)Ti 



Po{t) +pCn[po(T) +eo(r) 



£o(r)+p(C„ + CT)[po(r)+eo(r)] 



(C12) 
(C13) 

(C14) 
(C15) 
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from which we can write the equations for 5N'^ = V ^^^n^J^po and 5E^'^ = V ^^^n^^^eo 

d 



5N^' = -2pU5N''' -pCnSE 



M 



dT 
d_ 
dT 



'-pM 



6E^' = -p{Cn + Ct)SN''' - p(2Cn + Ct)SE 



M 



(C16) 
(C17) 



As we obtain a Hnear system of coupled equations, it is convenient to introduce some new variables to diagonalizc the 
problem: 



Cn 



6E 



M 



= + 6E 



M 



'-pM 



We obtain 



X. 



M 



-p(3C„ + Ct)X 



M 
2 • 



(C18) 
(C19) 

(C20) 
(C21) 



Let us suppose now that the equations for Xi{T) have the form of the macroscopic equations (|C20[) and (|C21[) plus 
a Gaussian random noise which variance has to be computed to reproduce the results that we have obtained for the 
correlation function of the fluctuations of N and E. The equations for the fluctuating variables Xi and X2 are thus 



d_ 



-pCn 



d_ 



X2 



i?2(r). 



(C22) 
(C23) 



Then, if we suppose that the noise terms have zero mean and its correlation function is delta-correlated in time 

{R^iT))H = 0, {R,iT')R,{T))H = ry5(T' - r), (C24) 

we can obtain the values of the amplitudes of the noise term in the same way we have done with the momentum, 
that is comparing with the results we have obtained in the previous section. We obtain 



(1 + z)'^ 
Til = -8pCJ—^An, 

r22 = -8p(3C« + Ct)(i + ^)'^22, 
ri2 = -4p(4C„ + Ct) ^ ' a,2. 



(C25) 
(C26) 
(C27) 



These calculations show that it is possible to describe the dynamics of fluctuations in the HDS in terms of some 
Langevin equations with Gaussian white noises. Due to the Gaussian nature of the noises and given that the equations 
are linear, the probability distribution function for those processes will be also Gaussian in agreement with our 
simulations. It is worth pointing out that, although the amplitude of the noises are known, they are not related in a 
simple way to the cooling rates C„ and (^t that appear in the "deterministic part" of the equations. 
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